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Abstract 

A discrete Laplace transform and its inversion formula are obtained by using a quadrature 
of the continuous Fourier transform which is given in terms of Hermite polynomials and 
its zeros. This approach yields a convergent discrete formula for the two-sided Laplace 
transform if the function to be transformed falls off rapidly to zero and satisfy certain 
conditions of integrability, achieving convergence also for singular functions. The inversion 
formula becomes a quadrature formula for the Bromwich integral. This procedure also 
yields a quadrature formula for the Mellin transform and its corresponding inversion 
formula that can be generalized straightforwardly for functions of several variables. 
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1. Introduction 



It is commonly accepted that the problem of obtaining a discrete formula for the 
Laplace transform 



of a function f(t) it is not so complicated as the inverse problem. This is due to the 
fact that the problem concerning the inversion of the Laplace transform is an ill-posed 
problem [5]. Many papers have been written on this subject (see [7j, which contains a 
list of references) but the techniques used to obtain an inversion formula can be classified 
PQ into only four main groups: those that use Fourier series, Laguerre functions, Gaver 
functionals, and the ones that discretize the Bromwich contour. The inversion formula 
for the two-sided Laplace transform presented here belongs to the last group and it is 
based on a quadrature of the integral Fourier transform |3j [2] . This quadrature formula 
is given in terms of a matrix N x N whose elements are constructed from the N zeros of 
the Hermite polynomial H^it) and has order 0(1/N) if the function to be transformed is 
square-integrable in (—00,00) and satisfy certain conditions of integrability [2j. The aim 
of this paper is to show that a simple and straightforward adaptation of such a formula 
yields a discrete two-sided Laplace transform with an easy-to-compute inversion formula 
corresponding to a quadrature of the Bromwich integral, and a discrete Mellin transform 
and its inversion formula. All of these discrete transforms can be generalized easily to the 
case of several variables. 



2. A Discrete Laplace transform 

Firstly, we reformulate the procedure followed in [3], [2] to obtain a quadrature formula 
for the integral Fourier transform yielding a discrete Fourier transform. Proofs and further 
applications can be found in these references. 

Let us consider the set of functions u n {i) = exp(—t 2 /2)H n (t), n = 0, 1, . . ., where H n (t) is 
the nth Hermite polynomial. This set is closed in L 2 (— 00, 00) [5] and their elements are 
related by the recurrence equation u n+ i(t) + 2nu n _i(t) = 2tu n (t), which can be written 
as the eigenvalue problem TU = W, —00 < t < 00, where T n k = 5 n+ i i k/2 + {n — l)5 n< k+i, 
n, k — 1, 2, . . ., and U is the vector whose nth entry is u n -i(t). 
The Fourier transform of u n (t), denoted by v n (u), is given by 



and satisfy the recurrence equation v n+ i(u>) — 2nv n -\(uj) = —2iu> v n (u>), which can be writ- 
ten in the matrix form W V = —iuY, —00 < u> < 00, where W n fc = 5 n+ i^/2 — (n— l)S n! k+i, 




(1) 
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n, k — 1, 2, . . ., and V is the vector whose nth entry is v n -i(u). By writing the recurrence 
equations as matrix equations we can consider the eigenproblems corresponding to the 
principal submatrices of order N of T and W to generate sequences of ^-dimensional 
vectors U y V converging to U and V respectively when N — > oo and in this way, to gene- 
rate approximations to the functions u n (t) and v n (cu). First let us note that the diagonal 
matrix § whose elements are given by Sj k = >j2 k ~ 1 (k — T)\5j k , generates a symmetric 
matrix S _1 TS and a skew-symmetric matrix § _1 WS whose principal submatrices of order 
N, denoted by T and W , have elements given by T nk = \/n/25 n +i,k + \/Jp — l)/25 n>fe+1 
and W nk = \Jn/25 n+ i t k — -\/{n — l)/28 n ^+i, respectively. Now let us consider the finite 
eigenproblems of T and W: 

TU k = t k U k , WV k = cu k V k , k =1,2,..., N. 

The above recurrence equations and the Christoffel-Darboux formula can be used to find 
the functional form of the eigenvectors, and also to show that the eigenvalues t k and u k 
are both zeros of H N (x). Thus, the nth entries of the eigenvectors U k and V k are given by 

U nk = <p n -i(t k ), V nk = (-i) n ~Vn-i(^fe), n = 1, . . . , N, (3) 

where 

. ^ ;(N-1)\2 N -™- 1 H m (x) 



Nm\ #jv-i(z) 

By construction, T and W approach T and W respectively when N — > oo. Therefore, in 
this limit, the nth elements of U k and V k approach u n (t k ) and v n {oj k ) respectively, up to 
a constant factor. Since v n {uj) is the Fourier transform of u n (t), the linear transformation 
F which yields the vector V nk , k = 1, . . . , N, when it is applied to U nk , k = 1, . . . , N, 
corresponds to a discretization of the Fourier transform. This transformation is determined 
by the matrices U and V, whose kth columns are just U k and V k respectively. Since F 
satisfies the relation V T = FU T between the transpose matrices V T and U T , we get 

F = V T U. (4) 

The elements of the unitary and symmetric matrix F 

f - 2N ' 1 ( N - 1 ) 1 yH)" (f)ff(uj) 



satisfy 



for bounded tj and u k . Here, At = tj + ± — tj = n /\/2N is the Riemann measure that yields 
the quadrature formula 



F kj = -^(-iy+ k e- u ^ + 0{1/N), 
v2ir 



e-^f{t)dt = / e- s ^f(t)dt = V^J^i-^Fkjfitj) + 0(1/N) 

-oo J —OO _2 



(6) 



3 



for the integral Fourier transform of f(t) evaluated at u k and for the two-sided Laplace 
transform of /(£) evaluated at s k = iu k . The order of this formula holds whenever f(t) 
satisfies certain conditions of integrability |2j. If furthermore/ (t) is a causal function 



h(t), t > 
0, t < 0, 



/(*) 

equation ([6]) becomes a discrete formula for the Laplace transform of hit) 

g(s k ) = / e- s ^h(t)dt = V L kj f(tj) + 0(l/N), (7) 
Jo j=1 

where s k = ioj k and 

L kj = V2^(-iy +k F kj . (8) 

The generalization of this discrete transform to several variables is straightforward. Let 
g^s 1 , s 2 , . . . , s n ) be the n-dimensional two-sided Laplace transform of fit 1 , t 2 , . . . , t n ), i.e., 



4 „2 



/oo 
e- s -7(*V 2 , . . . ,t n )dt l dt 2 ■ ■ ■ dt n , 
-oo 



where s = (s 1 , s 2 , . . . , s n ) and t = (t 1 , t 2 , . . . , t n ). Then, the corresponding discrete trans- 
form if given by the matrix 

L = L n ® - • • ® Li ® • • • ® Li (9) 

in which the entries of L\ are built out of Ni Hermite zeros lying on the Zth direction and 
the approximant g to ^(s 1 , s 2 , . . . , s n ) is obtained through the product 

g = Lf, (10) 

where L is the matrix defined in ([9]), f is the vector whose components are given and 
ordered by 

fr /('!,•'},•••• -'I)" (11) 

The index r is related to the others by r = j\ + (j 2 — l)Ni + (j 3 — l)NiN 2 + ■ ■ ■ + (j n — 
1) YTi=i Ni, where ji = 1, 2, . . . , iVj. The component g r of the vector g is the approximation 
to the exact transform g{s^ v s 2 2 ,--- , s"J where s*- = iu l j , = ^- y / = 1, 2 . . . , n. 
Notice that this approach on the discrete Laplace transform put the direct and inverse 
problem on the same footing since F~ l = F'. Thus, the elements of L~ l can be computed 
directly by 



L- k * = (-l) j+k F* k /V2 



7T, 



where * means complex conjugation. By applying L^ 1 to the vector whose elements are 
the values of g(s) evaluated at s k = iu k , we get an approximation to the values of fit) at 
tj. In other words, we obtain a discretization of the Bromwich integral 

-j r-ioo N 

fitj) = 77- / e^'g(s)ds = J2 L Jk9(s k ) + 0(1/N), s k = iu k , (12) 
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in which the contour of integration is the imaginary axis and the singularities of g(s) lie 
on the left of this line. The real part and the imaginary part of g(s) should satisfy the 
conditions on integrability given in [5]. The extension to several variables is obvious. Now 
we have 

f^L^g, (13) 

where 

L 1 = L^ 1 ® - • • ® L~[ l ® - • • (g) L^ 1 , (14) 

g is the function g^s 1 , s 2 , ■ ■ • , s n ) evaluated at s*- = iu l j , I = 1, 2...,n, and f is the 
approximant to /(i 1 , t 2 , . . . , t n ). 



3. A discrete Mellin transform 

Since the Mellin transform 

/■oo 

9m(s) = / x 8 ~ 1 f(x)dx 
Jo 

is a two-sided Laplace transform under the transformation x = exp(— t), the discrete 
Laplace transform L defined in ([9]) yields a discretization of the multidimensional Me- 
llin transform evaluated on the imaginary axis of each variable s l , I = 1, . . . ,n. Thus, 
we have that if fix 1 , x 2 , . . . , x n ) is the function to be transformed, the pair of discrete 
multidimensional Mellin transforms are given by the formulas 

g M = Lf t , f t = L^gM, (15) 

in which f t is the vector whose elements are given by 

f t (t\t 2 , ...,t n ) = /(exp(-t 1 ),exp(-t 2 ), . . . ,exp(-r)) 

and ordered according to (TTTT) . It should be noticed that in the inverse formula, the vector 
f t approaches ft(t k ) instead of f(x k ). 



4. Examples 

In this section we perform some numerical calculations to show the accuracy of the 
above discrete Laplace and Mellin transforms. We present two singular cases (the first 
and third examples) for which the discrete transforms yield convergent results. For such 
cases the necessary conditions to get the order 0(1/N) are not fulfilled, therefore, the 
order of convergence is estimated numerically in the next section. 
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4.1. Discrete Laplace transforms 



As a first example, we compute the numerical inversion of 

n 

g (s) = 2^cosh(fcs), (16) 

k=l 

which is the two-sided Laplace transform of a train of 2n delta functions centered at the 
integers ±1, ±2, . . . , ±n. This problem resembles the numerical inversion of the partition 
function of the quantum harmonic oscillator, a typical test problem. 
In order to approximate the inverse transform of ffTBl) . the number N of zeros of Hff(t) 
should be greater than n 2 /2 because in this way the interval [— n, n] is contained in 
(—\/2N + 1, \/2N + 1), which is the interval where the Hermite zeros lie. The application 
of L^ 1 to the vector g whose elements are the values of ( [TBI at the Hermite zeros on 
the imaginary axis yields the interpolated set of points shown in Figure 1. The result is 
a function showing the typical features of a sum of delta functions centered at integer 
values. 




-6 -4 -2 2 4 6 -10 -5 5 10 



♦—Real Part * — Imaginary Part 



Figure 1: Two-sided numerical inversion of the Laplace transformed function 
(|16D . In (a) 40 Hermite zeros have been used for n = 4. In (6) 100 Hermite 
zeros have been used for n = 13. The maxima of the real parts are centered 
at the corresponding integers and the imaginary parts are zero. 



As a second example, we take the function h(t) = exp(— t) sin(t), < t < oo, whose 
one-sided Laplace transform is given by g(s) = l/[(s + l) 2 + 1]. According to ©, h(t) 
should be substituted by the causal function 

f/ t \ = f ex PH) sin (*)> 1 ^ (17 \ 

u |0, t<0, 1 ' 



in order to obtain the approximated Laplace transform. The application of (JTj) and ( jl~2l) 
to the vectors / and g respectively, yields the results displayed in Figure 2. For iV = 40, 
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the relative errors are given by 

\\g-gh 



0,023758, 



11/ -/II: 



0,0236836. 



It should be reminded that g, g and / are complex vectors. 
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■Exact function 



■Approximate result 



Figure 2: (a) Exact and approximate Laplace transform of (|17f) . (b) Inver- 
se transform obtained through (|13p . In both cases, 40 Hermite zeros on the 
imaginary axis were used. 



The next examples concern the performance of the discrete Mellin transform. As a first 
case, we take the singular problem defined by the Mellin transform of 



fix) 



X 



1 — x' 



< x < oo. 



The Cauchy principal value of this integral is — 7rtan(7rs) and it is displayed in Figure 3, 
together with the discrete Mellin transforms (jl~5l) . Figure 3(6) shows the plot of /(exp(— £)) 
against t instead of f(x) against x. The corresponding relative errors are 



\9m — 9m\\2 

||<7iW 2 



0,156919, 



Wft- ft 

11/ 



t 2 



0,0739943. 



t 2 
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Figure 3: (a) Exact and approximate Mellin transform of (I18p on the imaginary 
axis (the real part is zero). (6) Inverse transform obtained by (|15p . In this case 
the imaginary part is zero. In both cases 40 Hermite zeros were used. 



As a final example, we consider the function 

f{x) = exp( 



< x < oo, 



(19) 



whose Mellin transform is sin(7rs/4)r(s). Figure 4 shows the output of the discrete trans- 
forms ( jl~5i) . Again, /(exp(— £)) is plotted against t in Figure 4(6). The relative errors are 



\9-9h 

\\gh 



0,00702041, 



11/ -/II: 



0,00701767. 
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Figure 4: ( a) Exact and approximate Mellin transform of (I19j) on the imaginary 
axis. (6) Inverse transform obtained by (I15p . In both cases 40 Hermite zeros 
were used. 
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5. Final remark 



Finally we address the performance of the discrete Laplace and Mellin transforms 
presented here on singular problems such as the above first and third examples. Repeated 
numerical calculations can be done to estimate the convergence of the results yielded 
by these discrete transforms. Thus, by changing the number N of Hermite zeros it can 
be seen that in the third example, the Mellin transform of (fl8"]l . the relative error goes 
as 1/yN. In the case of the first example, the Laplace inversion of f|T6|) . it is necessary 
to measure convergence in a different way since it is not possible to evaluate a delta 
function. To this end, we compute the area under the linear interpolation of the entries 
of the vector yielded by the numerical Laplace inversion, and test this value against the 
correct result. For n = 1 this integral should be 2 and the numerical integrations give 
2.0052, 2.0032 and 2.0025, for 50, 80 and 100 Hermite zeros, respectively. In order to give 
a visual representation of this case, we present in Figure 5 the discrete inverse for n = 1 
and N = 100. 




-10 -5 5 10 



Figure 5: Numerical inversion of the Laplace transformed function (|16p for 
n = 1 and N = 100. 
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